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Grand canonical simulations designed to resolve critical universality classes are reported for z:\ 
hard-core electrolyte models with diameter ratios \ — a+/a-<6. For z = l Ising-type behavior 
prevails. Unbiased estimates of T C (X) are within 1% of previous (biased) estimates but the critical 
densities are ~5 % lower. Ising character is also established for the 2:1 and 3:1 equisized models, 
along with critical amplitudes and improved T c estimates. For z = 3, however, strong finite-size 
effects reduce the confidence level although classical and 0(n> 3) criticality are excluded. 
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Since the challenging experiments on Coulomb- 
dominated electrolyte solutions more than a decade ago, 
the nature of ionic criticality has been an outstanding ex- 
perimental and theoretical issue Q. A central question 
has been: Is ionic criticality of Ising-type as for simple, 
nonionic fluids, or of mean-field character, reflecting the 
long-range Coulombic interactions, or something still dif- 
ferent? To tackle this problem theoretically, hard-sphere 
ionic fluids (the so-called primitive models) have been 
under broad study analytically and via simulations. 

These models consist of N = N + +N- hard sphere ions 
in a volume V , N + of diameter a + and charge q+ = zq 
and N- of diameter a_ with charge q_ = — q. Ions a and 
r interact via the Coulomb potential q a q T /r. Extensive 
simulations of the simplest model, namely, the restricted 
primitive model (RPM) with z — 1 and a + = a_, together 
with appropriate finite-size scaling analyses, have estab- 
lished unequivocally that the critical behavior belongs to 
the expected (d = 3) Ising universality class with criti- 
cal exponents /3~0.326, 7 ~ 1.239, etc. HHH. How- 
ever, the RPM possesses artificial symmetries of both 
charge and size. Thus, to better understand ionic criti- 
cality, an important question has remained open. Specif- 
ically, does the universality class remain unchanged when 
the symmetry of size or charge is broken? Although 
experiments on real electrolytes (which are, of course, 
size-asymmetric) favor Ising-type criticality 0, the 
presence of short-range forces between ions and solvent 
molecules tends to cloud the theoretical issues regarding 
the effects of size and charge asymmetry on ionic criti- 
cality. 

Now long-range Coulombic interactions are known to 
be screened exponentially at high temperatures and low 
densities. But whether exponential screening prevails at 
criticality is still an open question. How, indeed, does the 
diverging density correlation length influence charge cor- 
relations near criticality? It has been held that the RPM 
maintains full screening even at criticality |(J, thereby 
supporting expectations of Ising character. On the other 
hand, Stell |6( has argued that when size or charge sym- 
metry is broken, the density fluctuations play a crucial 



role, mix into the charge correlations, and cause the 
charge screening length to diverge on approach to criti- 
cality. He further suggested that this might change the 
critical universality class of an ionic fluid leading to mean 
field behavior 0. 

Aqua and Fisher Q addressed this issue by analyzing 
exactly soluble ionic spherical models. For z = l, they 
proved that, indeed, when asymmetry is present (which 
can arise from short-range interactions yielding, in effect, 
a + 7^a_) the density-density and charge-charge correla- 
tions mix and the charge correlation length diverges at 
criticality, as concluded by Stell Q- Nevertheless, con- 
trary to Stell's proposal |g, the spherical-model critical- 
ity is preserved. The coexistence curves of spherical mod- 
els are, however, always parabolic (i.e., (3 — i) in contrast 
to the RPM. Furthermore, spherical models with z > 1 
have not as yet been studied. Hence it remains a chal- 
lenge to resolve the question of the critical behavior of 
asymmetric primitive models. 

Some years ago simulations of nonsymmetric primi- 
tive models were performed H, 0, 0] . It was found that 
suitably normalized critical temperatures (see below) de- 
crease with size and charge asymmetry while the critical 
densities increase with z but decrease with size asymme- 
try. To compute these critical parameters, however, the 
Bruce- Wilding (BW) method [ll| was employed: this as- 
sumes from the outset the Ising character of criticality. 
Thus, although the simulations were consistent with Ising 
behavior, they did not rule out other possibilities, such 
as classical, i.e., van der Waals (vdW), or selfavoiding 
walk (SAW), or XY (i.e., rt = 2), etc. Furthermore, no 
allowance was made in the analysis for pressure mixing 
in the scaling fields [l2j . 

Our aim here, therefore, is to determine from first prin- 
ciples the universality class of asymmetric primitive mod- 
els and, further, to obtain unbiased estimates for the crit- 
ical parameters taking adva ntag e of recently developed 
finite-size scaling techniques [l3| . 

Accordingly, we have performed grand canonical 
Monte Carlo (GCMC) simulations on z:\ primitive mod- 
els. Reduced temperatures and densities are defined 
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naturally by T* = k B Ta±/zq 2 and p*=Na\/V where 
a± = ^(a+ + a_) is the unlike-ion collision diameter [Io| . 
The size asymmetry is described by A = a + /a_ • For 
size asymmetry (with z~l), we compare results for A = 3 
and 5| with the RPM (A = 1); for charge asymmetry we 
study z = 2 and 3 with A = 1. 

Our calculations have been eased by considering 
finely-discretized primitive models at discretization level 
(=a±/ao = 10 so that particles are allowed to sit only on 
lattice sites of spacing ao 0] . Continuum models corre- 
spond to the C = °o limit. However, it has been shown 
that choosing £ > 4 does not change the universal critical 
behavior and makes only small changes in the critical pa- 
rameters 0| . Thus there is little doubt that criticality 
in (C = 10)-level models belongs to the same class as in 
the continuum. To maximize the information gleaned, 
multihistogram reweighting techniques have been used. 

Following our previous strategies for the RPM 0,0, 0' 
we have studied the finite-size moment ratios 
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where (-)l denotes a grand-canonical expectation value 
in a periodic box of side L at fixed T and chemical po- 
tential p chosen to yield the desired mean density (p) L ■ 
Now, in a single-phase region Q L — > | when L — > oo, while 
Ql — ► 1 on the coexistence-curve diameter. At criticality, 
Ql(T c , Pc), approaches a universal value, Q c , that serves 
to resolve distinct universality classes of criticality rather 
sharply. Thus, for Ising systems one has Q\ 8 — 0.6236, 
while the vdW, SAW and XY values are 0.4569 • • • , 0, and 
0.804 5 , respectively; for long-range, l/r 3+cr Ising systems 
with | < a < 1.966, Qc{<*) rises smoothly from 0.4569 • • ■ , 
reaching Q c ~ 0.600 at a= 1.9 0. 

To go further it proves important to calculate Q L along 
the Q-loci, pq(T;L), in the (p,T) plane defined by the 
values of (p) L for which Ql is maximal at fixed T [2l.ll3|: 
at T = T C these loci approach p c when L—>oo. Plots 
of Ql(T) evaluated on the Q-loci for sufficiently large 
'nearby' sizes, L— 1 and L, cross one another near T c , say 
at T®(L), with a value, say Q®(L), close to the corre- 
sponding universal value Q c . Finite-size scaling theory 
[l3j implies that Q®(L) approaches Q c as L~ e l v ', while 
Tf?(L) approaches T c more rapidly, as 1/ L^ 1+e ^ v . For 
Ising universality, one has i/~0.63 and #~0.52. 

Data for size-asymmetric primitive models with 
A = 3 and 5| are shown in Fig. ^ for sizes from 
L* = L/a± = 6 to 13 [lfj. As for the RPM (see Fig. 1 of 
Ref. 0), the self-intersection points almost coincide with 
the universal Ising value, Q l c s . Furthermore, Q^ dw and 
Q^ Y are off scale while a limit at Q c < 0.600 is implau- 
sible. Thus these plots not only confirm Ising-type be- 
havior but also exclude vdW, SAW, XY, and long-range 
Ising criticality with a < 1.9. Notice, however, that the 
convergence of Q®(L) becomes slower when A increases; 
so a convincing study for larger values of A would require 
still bigger systems. Nevertheless, we conclude that size 
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FIG. 1: Plots of Ql on the Q-loci vs T* for size asymmetric 
1:1 primitive models. The horizontal line marks Ql = Qc S ; 
the arrows locate estimates of T* . 



asymmetry has no effect on the character of ionic criti- 
cality in accord with the spherical model results 0). 

On the other hand, all nonuniversal critical parameters 
depend on A Previous estimates ofTU and 

p c were obtained via the biased BW procedure [ill Il2|: 
but the plots of Fig. ^ yield the unbiased precise values 
presented in Table ||| These values for T* are only 0.2% 
lower than the earlier estimates except for a surprising 
~ 1% discrepancy when A = 3; but our results are more 
precise by factors of 5 or more. 



TABLE I: Monte Carlo estimates for the size asymmetric 1:1 
primitive models at discretization level £ = 10 |4. ll4lHq| . The 
uncertainties in parentheses refer to the last decimal place. 
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4.952(2) 7.6(2) 
4.472(2) 6.1(2) 
3.654(2) 4.3(2) 


4.96(1) 7.9(2) 
4.42(1) 6.4(4) 
3.66(8) 4.6(1) 



To obtain the estimates for p* listed in Table [I] we 
have, as previously 0, extrapolated the Q-loci values, 
p Q {T*;L) vs l/Lf 1 -")/" with a ~ 0.11. The original BW- 
based estimates are higher than ours by 5, 4 and 7% (in 
accord with RPM studies 0). These discrepancies may 
be due to significant pressure mixingll2j. Overall, we 
confirm the previously discovered 00] strong decrease 
of T*(X) and p*(A) when A increases. 

The charge-asymmetric primitive models are of 
especial interest and pose a greater challenge owing to 
strong ion-association in the critical region leading to a 
jopulation of tightly bound neutral and charged clusters 
lS l 1 191 . The models have been studied by simulation 
9j, UJJ, yjj and recently an effective Debye-Hiickel-type 
theory has been developed 0. One finds that T*(z) de- 
creases steadily when z increases while p*(z) rises sharply 
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liilliiiLLi! ■ Although there is a general belief (adopted 
in the simulations [ii UJJ, llH ) that increasing z does not 
alter the Ising critical behavior, there has been no con- 
vincing supporting evidence. 

Our simulations of the equisize 2:1 and 3:1 hard-core 
models address this issue jlfj|. The successive intersec- 
tions, Q$(L) and T C Q (L), for z = 2, found as in Fig. □ 
are presented in Fig. [3 for L* > 11. In contrast to the 
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FIG. 2: Plots of Q?(L) and T C Q (L) vs (L* -lo)^ with V = 1 
and (1 + 9) jv ~2.41, respectively, for the equisized 2:1 elec- 
trolyte. The shifts lo — 0, 1, 2, and 3 (plotted with circles, 
diamonds, triangles and crosses) are merely aids to extrapo- 
lation: the dashed lines are guides to the eye. 

RPM (see Fig. [T] and Refs. H0), the self-intersections 
converge much more slowly to their limits as L increases. 
Nevertheless, the plots in Fig. j2{a) strongly favor Ising 
criticality and surely rule out the 'adjacent' vdW and 
XY possibilities at Q c ~ 0.457 and 0.805. We also find 
T*(z = 2) = 0.0466(1) and p* c (z = 2) = 0.093(3); the for- 
mer is 0.9% lower than the estimate 0.047, obtained by 
using the BW method at L* = 15 10], while the value for 
p* agrees [io| . 

To understand the behavior of (L) in further detail, 
we also simulated smaller systems, L* < 12, for z = 2 and 



3 20] . Surprisingly, the behavior of the Q-intersections 
turns out to be sharply nonmonotonic in L: see Fig- EI In 
fact, this is also seen for the RPM (z = 1) when data for 
L* < 6 are examined. Although unexpected, the behav- 
ior is clearly a reflection of strong finite-size effects. To 
gain a feeling for the phenomenon, note that in GCMC 
simulations, one attempts to insert a neutral cluster or 
"molecule" of z+1 ions. When z increases, larger systems 
are evidently needed to accommodate the same number 
of molecules at comparable molecular densities. One may 
thus enquire into the geometry of neutral ion clusters. To 
reflect this, the data in Fig. |3] have been plotted using a 
system size, L, rescaled by za± (open symbols) and by 
y/z + la± (solid symbols). The first scaling embodies 
the center-to-center size of a neutral cluster on a line; 
indeed, for z = 1 and 2 the molecular ground states are 
linear. But for z = 3 the ground state is a planar centered 
triangle 0: the second rescaling respects this feature 
and seems more appropriate. 
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FIG. 3: Plots of Qg(L ) vs with V = L* / z (open 

symbols) and L* / 'V ' z + 1 (solid symbols) for z = l, 2 and 3. 

Also relevant and seen in Fig.|3|is the sharp changeover 
in behavior with a strong z-dependence at relatively 
small sizes. Finite-size scaling theory ^^] predicts 
{Ql{T®{L)) - Q c ] - L- Q l v with higher order terms 
L~ 2/3 / u , L~( 1 ~ Q 0/ l/ , etc. The interference of these terms 
may certainly lead to nonmontonic behavior; however, 
naive fits to the data with a few leading terms are not 
successful. Furthermore, with only sizes up to L* = 11 
available one would be led to guess that the limiting Q c 
for z = 2 and 3 exceeded 0.8 or 0.9, respectively; but, as 
clear from the z = 2 data, this is not sustainable. Nev- 
ertheless, on the basis of Fig. |3 it is tempting to spec- 
ulate that the z > 2 equisized primitive models may ex- 
hibit some sort of crossover from n = oo behavior (with 
~ 1.0) for z ^> 1 to n = 1, Ising behavior for finite z as 
L increases. But whether this is a genuine effect surely 
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needs further investigation! 

In the absence of precise data for significantly larger 
systems [2(J, we cannot convincingly argue that the 
3:1 equisized model exhibits Ising behavior. Thus, 
while vdW criticality seems excluded, the turnover in 
Fig. El for L* = 12, is consistent with XY or 0(n = 2) 
criticality. Nevertheless, by comparison with the 
z = 2 data (which extend to L* — 16), our belief that 
Ising criticality will still be realized when L* — > oo 
is surely reasonable. Accepting that, one may esti- 
mate T* by extrapolating the intercepts of Ql(T*) 
with Qjf; then p*(z — 3) follows from the Q-loci, as 
for z=l and 2. This yields T*(z = 3) =0.0405(2) 
and p*(z — 3) = 0.126(3) which may be compared with 
0.0410(1) and 0.125(4) J respectively, obtained by the BW 
procedure at L* = 15 [lOj . 

Finally, the critical amplitudes are of interest j2]| . For 
the RPM we have C + a\~ 0.087 for the susceptibility 
and Ba± ~ 0.284 for the order parameter j^]. Together 
with the universal ratio Q c = (tf, ) d B 2 /C + ~ 0.323 
[2ll | these yield the density correlation-length ampli- 
tude, £+ /a±~ 0.704; this is surprisingly close to 
the generalized-Debye-Huckel (GDH) value 0.7329 [23| . 
For z — 2 and 3, we now find, to lower precision, 
C+a\ = 0.11(1), 0.16(3) and Ba\ = 0.35(5), 0.45(5) £3, 
respectively, which yield £i Q /a± = 0.66(5) and 0.63(9), 
not inconsistent with the RPM value. The relative in- 
sensitivity of ^ to z is notable. 

In summary, we have studied the critical behavior of 
size- and charge- asymmetric primitive model electrolytes 
via extensive Monte Carlo simulations. It has been 
shown convincingly for z < 2 that breaking either symme- 
try leaves the universal Ising character unchanged. For 
the equisized charge-asymmetric 3:1 electrolyte, strong 
finite-size effects mandate simulations of much larger sys- 
tems; nevertheless, the data are consistent with Ising- 
type behavior and van der Waals and 0(n > 3) criticality 
are excluded. The density correlation-length amplitudes 
vary little with z and are to be close to the prediction of 
GDH theory 0. 
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